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ABSTRACT 

Supermassive black holes with up to a 10 9 M© dwell in the centers of present-day galaxies, 
and their presence has been confirmed at z > 6. Their formation at such early epochs is still 
an enigma. Different pathways have been suggested to assemble supermassive black holes in 
the first billion years after the Big Bang. Direct collapse has emerged as a highly plausible 
scenario to form black holes as it provides seed masses of 10 5 - 10 6 M . Gravitational collapse 
in atomic cooling haloes with virial temperatures T vir > 10 4 K may lead to the formation of 
massive seed black holes in the presence of an intense background UV flux. Turbulence plays 
a central role in regulating accretion and transporting angular momentum. We present here 
the highest resolution cosmological large-eddy simulations to date which track the evolution 
of high-density regions on scales of 0.25 AU beyond the formation of the first peak, and 
study the impact of subgrid- scale turbulence. The peak density reached in these simulations 
is 1.2 x 10" 8 g cm -3 . Our findings show that while fragmentation occasionally occurs, it 
does not prevent the growth of a central massive object resulting from turbulent accretion and 
occasional mergers. The central object reaches ~ 1000 M© within 4 free-fall times, and we 
expect further growth up to 10 6 M© through accretion in about 1 million years. The direct 
collapse model thus provides a viable pathway of forming high-mass black holes at early 
cosmic times. 
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1 INTRODUCTION 



The formation of pregalactic black holes is a fundamental issue 
| in the study of galaxy formation. We know that most galaxies 
have black holes in their nuclei. Numerous quasars with black hole 
masses of 10 9 M© have been observed at redshifts higher than six 
when the Universe was less than a billion years old (IFan et al.l2003l 
2006). Recently, the discovery of a z=7.1 q uasar with a black hol e 
mass of 2 x 10 9 M is even more striking (IMortlock et al.ll201ll) . 
It is however not very well understood how these quasars observed 
in the Sloan Digital Sky Survey could attain such high masses in a 
such short span of time after the Big Bang. It suggests that black 
hole seeds have been formed even earlier and it is important to get 
a physical understanding of their origin. 

One scenario for the formation of supermassive black holes 
(SMBHs) could be that th ey formed throug h accretion and merg- 
ing of stellar remnants (JHaimanl l2004bU ah. Numerical studies 
show low gas densities in the vicinity of black holes remnants 
due to their feedback which brings an obstruction for their for- 
mation though accr etion and merging (jjohnson & Bromml 120071 ; 
lAlvarez et al.l 12008). Black holes may have formed through run- 
away collisions of young dense stel l ar clu s ters due to re l ativis- 
ts instability ( Portegie s Zwart et al.l 120041: bmukai et all 120081 : 
iDevecchi & Volonte ri 2009) but this can happen only if the Uni- 
verse is enriched with trace amounts of metals after the for- 



mation of the first stars between redshift 10-20, and the result- 
ing black hole masses are about 10 3 M©. Primordial black holes 
may have formed during the early stages of the Big Bang al- 
beit we do not see any observational evidence of their formation 
flRicotti et al.ll2008r) . Various other w ays for the formation of mas- 
sive b l ack holes has been s ugges te d dRee s 1978; Diorgovski et al.l 
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Johnson et alj2013l : lReisswig et all2013b . Protogalactic metal free 
haloes may directly collapse to form supermassive black holes pro- 
vided that the fragmentation in these halo s remains suppressed 
dOh & Haim an 2002; B romm & Loebll2003l: ISoaans & SfflJ |2006: 
Begelma net alJl2006l: iDiikstra et al 1120081: iDiorgovski et alJ l2008: 
Sha ng et alJbOld : ISchleicher et al.ll2010l : lLatif et aLlboTTIT Given 



the difficulties with other methods for the formation of supermas- 
sive black holes, direct collapse emerges as a plausible scenario. 
Haloes with virial temperatures > 10 4 K are the most likely can- 
didates for the formation of such objects. The formation of black 
holes via direct collapse is feasible if fragmentation is sufficiently 
suppressed, so that most of the material is accreted onto the central 
object. 

The essential condition for the direct collapse scenarios is to 
avoid fragmentation. Molecular hydrogen is the only coolant in 
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metal free haloes which can bring the temperature of the gas down 
to a few hundred Kelvins and may induce fragmentation by reduc- 
ing the Jeans mass. The gas is unlikely to fragme nt to form star s 
in primordial halos in the absence of H 2 cooling dLi et al.ll2003l) . 
For zero metallicity atomic cooling haloes devoid of molecules, the 
Lyman alpha radiation is the only pathway to cool the gas down 
to 8000 K. Cosmological simulations performed to study the col- 
lapse of protogalactic metal free halos show that self gravitating 
gas collapses almost i sothermally in the absence of H 2 coolin 



(iBromm & Loebl2003l : IWise et al.l2008l : lRegan & HaehnelJ2009bl : 
lLatif etalJEoTTh . 

In this respect,the presence of UV radiation has important 
implications for the formation of supermassive black holes as 
it may suppress cooling and consequently halt the fragmenta- 
tion of a gas cloud. It is found that strong Lyman Werner flux 
> 10 3 in units of J 2 i is required to photo-dissociate H 2 molecules 
in haloes with virial temperatures > 10 4 K ( fOmukail 12001 
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Agarwal et al1l2012l : lLatif et al.ll2013al) 

Turbulent energy plays a key role in the formation of struc- 
tures as it locally compresses the gas to high densities where it 
may become susceptible to fragmentation and collapse. It can pro- 
vide an additional pressure against the gravitational collapse at 
larger scales. In order to transport the angular momentum, previous 
studies suggested sce narios based on gravitati onal instabilities in 
self-gravitating disks (JKoushiappas et al. 2004 J) and self -regulatio n 
due to the presence of turbulence ( Begelman & Shlosm an 2009). 



In a similar fashion, turbulence is th ought to regulate the angu- 
lar m omentum transport in minihalos (I Abel et al.ll2002l; iTurk et all 
120121) . In the presence of high mass accretion rates and efficient 
cooling, gravitational torques may not be sufficient to transport 
the gas through the d isk leading to fragmentation (JClark et al.l 
1201 1| ; iGreif et al.ll2012l) . Numerical simulations show that primor- 
dial gas becomes turbulent during th e formation of the first galaxies 
dWise et al.ll2008l : IGreif et~aDl2QQ8h . While these studies employed 
a typical Jeans length resolution of 16 cells, later studies demon- 
strated the necessity of a higher Jeans resolution i.e., minimum res- 
olution of 32 cells per J eans length to obtain converged tur bulen t 
energy (iFederrath et al.11201 lb ITurk et al.lEoi2L lLatif et al.lEoi3ah . 
In addition, the turbulent energy cascade is expected to generate 
turbulence even on subgrid- scales (SGS), which may backreact on 
the resolved scales via turbulent viscosity. lLatif et al.l (I2013ah have 
shown that applying a SGS turbulence model to cosmological AMR 
simulations significantly changes the morphology of the halo dur- 
ing its initial collapse phase in comparison to standard simulations 
without an explicit SGS model. This may have a significant impact 
on the the formation of astrophysical objects. 

The main objective of this work is to assess the role of turbu- 
lence in the formation of seed black holes. To accomplish this goal, 
we have performed the highest resolution cosmological large eddy 
simulations to date which track the gravitational collapse beyond 
the formation of the first peak in the presence of a strong Lyman 
Werner background UV flux. We employ a fixed Jeans length of 64 
cells throughout the evolution of simulations to resolve turbulent 
eddies and follow the collapse down to sub AU scales for 4 free-fall 
times by exploiting the adaptive mesh technique. This study will al- 
low us to explore the feasibility of the direct collapse scenario. It 
will further enable us to quantify the impact of subgrid/resolved 
turbulence in the formation of supermassive black holes. 

Our paper is organized as follows. In the next section, we de- 
scribe the simulation setup and the numerical methods employed. 



In the 3rd section of the paper, we present the results obtained in 
this study. In the last section of this article, we summarize our main 
results and confer our conclusions. 



2 BASIC SETUP 

2.1 Computational methods 

We present high resolution cosmological simulations which capture 
a large range of spatial scales starting from Mpc down to sub AU 
scales. These simulations were performed usin g a modified version 
of the adaptive mesh refinement code Enzo JO'Shea et al.l 12004) 
including a subgrid-sc ale (SGS) model wh ich takes into account 
unresolved turbulence dSchmidt et al.ll2006l) . Enzo is an adaptive 
mesh refinemen t (AMR), parallel, grid-based cosmologi cal hydro- 
dynamics code JO' Shea et al.ll2004l : iNorman et al.ll2007l) . It makes 
use of the message passing interface (MPI) to achieve portability 
and scalability on parallel systems. The computational domain is 
discretized using a nested grid structure. We use a split hydro solver 
with a 3rd order piece- wise parabolic (PPM) method for hydro dy- 
namical calculations. The dark matter N-body dynamics is solved 
using the particle-mesh technique. A multigrid Poisson solver is 
employed for the self- gravity computations. 

The simulations are commenced with cosmological initial 
conditions generated from Gaussian random fields. We make use 
of the inits package available with Enzo to create nested grid ini- 
tial conditions. Our simulations start at redshift z = 100 with a top 
grid resolution of 128 3 cells and we select a massive halo at red- 
shift 15 us ing the halo finde r based on a standard friends of friends 
algorithm JTurket al.ll201ll) . Our computational volume has a cos- 
mological size of 1 Mpc h _1 and is centered on the most massive 
halo. Two additional nested levels of refinement are subsequently 
employed each with a resolution of 128 3 cells. In all, we initial- 
ize 5767168 particles to compute the evolution of the dark matter 
dynamics yielding a dark matter resolution o f 620 M Q . The pa- 
rameters from the WMAP seven years data djarosik et al.l 1201 lb 
are used for creating the initial conditions. We further allow ad- 
ditional 27 levels of refinement in the central 62 kpc region of 
the halo during the course of simulation yielding a physical res- 
olution of sub AU scales (3 AU in comoving units). The reso- 
lution criteria used in these simulations are based on the Jeans 
length, the gas over-density and the particle mass resolution. The 
gas overdensity criterion flags grid cells for refinement when the 
gas density exceeds four times the mean baryonic density. Simi- 
larly, grids are marked for refinement when the dark matter den- 
sity exceeds 0.0625 times poMr la where r=2 is the refinement fac- 
tor, 1 is the refinement level and a = -0.3 makes the refinement 
super-Lagrangian. The grid cells matching these requirement are 
marked for refinement. We mandated a Jeans length resolution 
of 64 cells throughout the evolution. This criterion was applied 
during the course of the simulations to ensure a high resolution 
in the central c ollapsing region, sim ultaneously ensuring the Tru- 
elove criterion JTruelove et al.ll 19971) . When the highest refinement 
level is reached, the thermal evolution becomes adiabatic. This 
approach enables us to follow the evolution of structures beyond 
the formation of the first peak until they reach a peak density of 

1.2 x 10~ 8 g/cm 3 which corresponds to four free-fall times. Our 
study consists of 9 distinct halos selected from various Gaussian 
random seeds. The results from the SGS turbulence model are com- 
pared with normal runs. In total, we perform 18 simulations. 
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2.2 Chemistry 

According to Big Bang nucleosynthesis, primordial gas is mainly 
composed of hydrogen and helium. For zero metallicity gas de- 
void of molecules, atomic line cooling is the only pathway to cool 
the gas down to 8000 K. While molecular hydrogen can poten- 
tially cool below a few thousand K, the presence of strong UV 
flux emitted by the first gener ation of stars h a s important implica- 
tions for structure formati on (JOmukail 12001 1: lOmukai et al.l [20 08 ; 
IWolcott-Green et al.ll201lh . A photon with energy between 11.2 
and 13.6 eV is absorbed in the Lyman Werner bands of molecu- 
lar hydrogen and brings it into an excited state. The H 2 molecule 
then decays to the vibrational continuum of a ground state and is 
dissociated. This two step process is known as the Solomon process 
(I Abel et al]|l997h . The photo-dissociation reaction is given by 



H 2 +y^H* 2 



H + H. 



(1) 



The critical value of the flux can be obtained by balancing the H 2 
formation and dissociation time scales. The value of the critical UV 
intensity is higher for the larger halos as compared to the miniha- 
los and varies from halo to halo. It is found that a strong Lyman 
Werner flux > 10 3 in units of J 2 i is required to photo-dissociate 
H 2 molecules in halos wit h virial temperatures > 10 4 K for a radia- 
tion temperature of 10 5 K JOmukail200ll:IJorinson & Bromml2007l: 
Diikstra et al.l 120081: [Schleicher et al.l l2010l: IWolcott-Green et all 
201ll : lLatifetal.ll201ll : ISafranek-Shraderetal.ll2012l) . J.t is the 



background UV flux below the Lyman limit (i.e. 13.6 eV) in units 
of 10~ 21 erg cm -2 s _1 Hz -1 sr _1 . 

To study the thermal evolution of the gas, it is crucial to 
model the chemistry of primordial gas in detail. In order to in- 
clude the primordial non-equilibrium chemistry, the rate equations 
of H, H + ,He, He + , He ++ , e~, H~, H 2 , HJ are self-consistently solved 
in the cosmological simulations. We use the H 2 photo-dissociating 
background UV flux implemented in the Enzo code, and employ an 
external UV field of constant strength 10 3 in units of J 2 i. We pre- 
sume that such flux is generated from a nearby star forming halo 
jDiikstra et alj 12008) and is emitted by Pop III stars with a ther- 
mal spectrum of 10 5 K. We include all relevant cooling and heating 
mechanisms like collisional ionization cooling, radiative recombi- 
nation cooling, collisional excitation cooling, H 2 cooling as well 
as H 2 formation heati ng. The chemistry solv er used in this work is 
a m odified version of I Abel et al.l dl997h and lAnninos etal.1 (J1997J) 
(see lTurk et all J2012h for details). 

We note that for an optical depth r > 10 7 , the escape time 
of Lyman a lpha photons become s longer than the gas free fall 
time. In fact JSpaans & Silkl d2QQ6h suggested that at columns above 
10 22 cm -2 Lyman alpha photons get trapped inside the halo which 
suppresses the cooling and may lead to an adiabatic colla pse. The 
impac t of Lyman alpha t r apping was explored in detail bv lOmukail 
J200ll) : ISchleicher eTaD bold) , finding that cooling due to addi- 
tional processes becomes relevant in these cases particularly the 
two-photon decay (2s- Is transition) and the 3p-2s transition of 
atomic hydrogen. Effectively, the temperature evolution is then 
very close to the evolution obtained from optically thin Lyman a 
cooling. We also note that for a radiation temperature of 10 5 K, H 2 
is mainly dissociated by the Solomon process while f or the stellar 
spectrum of 10 4 K the main dis sociation route is H~ JShang et al.l 
l2010l : IWblcott-Green et al.l201ll) . In principle, one would of course 
expect the presence of both contributions. As we focus on a situ- 
ation where H 2 cooling is not relevant, the details of these pro- 
cesses are however not important here. We neglect here the effect 
of H 2 self- shielding which could potentially raise the strength of 



the Lyman Werner flux required to photo-dissociate molecular hy- 
drogen. The presence of X-rays/cosmic rays flux may boost the 
form ation of H 2 and further increase the th reshold of UV flux 
(ITanaka & Haim an 2009; Inavoshi & OmukafcOllI) . 



2.3 Turbulence Subgrid Scale Model 

The phenomenon of turbulence describes a process where energy 
cascades from large scales to small scales through a series of ed- 
dies. The formation and evolution of cosmic structures is a typical 
case of turbulence generation. Due to the high Reynolds numbers 
relevant for astrophysical systems, it is not possible to resolve all 
scales down to the dissipative scale even with adaptive mesh refine- 
ment techniques. In AMR simulations, mainly the core of the halo 
is well resolved and turbulence cascades from coarser grids corre- 
sponding to large scales down to the center of the halos without be- 
ing properly accounted for. To compute the unresolved turbulence 
on grid scales in our simulation, we perform large ed dy simulations 
(LES) with the SGS turbulence model proposed in ISchmidt et al.l 
(120060 . This SGS model is based on a mathematically rigorous ap- 
proach separating the resolved and unresolved scales, and connect- 
ing them via an eddy-viscosity closure for the non-linear energy 
transfer across the grid scale. The turbulent viscosity is given by the 
grid scale and the SGS turbulence energy, i. e., the kinetic energy 
associated with numerically unresolved turbulent velocity fluctua- 
tions. In contrast, implicit large eddy simulations (ILES) use only 
the numerical dissipation stemming from the discretization errors 
of the compressible fluid dynamics equations. This is the standard 
method in computational astrophysics. 

Applying the filtering mechanism to the fluid equations 
and solving them in a comoving coo rdinate system, we obtain 



ana solving mem in a comovmg coorc 
dSchmidt et alJ200d : iMaier et al J2009h : 



did 

T-P+ ~—VjP = 

dt a oxj 



(2) 



d 1 d 

Li- kJJvi 



dt 



d 1 d 

P^res + -T—Vjpen 

a oxj 



dt 



a dxj J a dxj lJ a J 



a oxi a a 3 



+ P(A + 6) - l a Vi £j Tij 



did 

T-pe t + -t— Vjpe t = . 
dt a dxj 



-p(A + e)- 



-Tij—Vt-2-pet 
a dxj a 



(4) 



(5) 



Here a(f) is the scale factor of the cosmological expansion, 
p and P are the comoving density and pressure, v t is the peculiar 
velocity, g* is the gravitational acceleration, £ res = e int + \v 2 the 
sum of the numerically resolve d thermal and k i netic e nergies, and 
e t the SGS turbulent energy fsee lSchmidt et al.l(l2006l) : lMaier et"aD 
(2009) for details). The terms D, A, e, and T tj correspond to SGS 
transport (diffusion), pressure dilatation, dissipation, and turbulent 
stresses. The following closure relations in terms of the resolved 
flow and the SGS turbulent velocity q = y2^t are applied to evalu- 
ate these terms: 



d ^ r 2 d 

= t-Cdp/a^ -r-q, 
dr t OYi 

r 2 9 
= C A q — v,-, 
dr t 



(6) 

(7) 
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= c f 



/a' 



\SijPq 



(8) 

(9) 

(10) 



where //» = pC v l A q is the dynamic turbulent viscosity and 

is the trace-free rate-of- strain tensor. Further details and numeri- 
cal validations of the above closures ar e given in dedicated studies 
dSchmidt et alfcOOd iMaier et alJ2009h . 

For our SGS model, the coefficients are calibrated against 
transonic compressible tur bulence simula t ions, c omparable to the 
regime in our simulations ISchmidt et al.l (120061) . The method of 
adaptively refined large eddy simulations is used to apply the SGS 
model in cosmological AMR simulations jMaier et al]|2009h . As 
the SGS turbulence energy depends on the grid scale, which varies 
on adaptive meshes, energy must be exchanged with the numeri- 
cally resolved velocity field when the grid is refined or de-refined. 
This is achieved by assuming the Kolmogorov two-thirds law for 
the scaling of the turbulent velocity fluctuations. As long as the tur- 
bulence is subsonic and nearly isotropic local ly, this is a reasonable 
assum ption. In contrast, the method used in iGrav & Scannapiecol 
(1201 lh does not calculate the turbulent energy on the grid scale, but 
the energy associated with a characteristic length scale of buoy- 
ant bubbles. Our SGS model completely neglects gravity on unre- 
solved length scales and assumes the turbulent cascade from larger 
resolved scales as the dominant source of unresolved turbulenc e 
for the effects of buoyancy on subgrid scales (ISchmidt et al.l2006h . 
With a resolution of 64 cells per Jeans length, we expect this as- 
sumption to be valid on unresolved scales. 

Even though the presence of turbulence was r eported in mas- 
sive prim ordial halos already bylwise et al.l (1200 8), it was only re- 
alized by iFederrath et al.l (J201ll) that a minimum resolution of 32 
cells per Jeans length is required to obtain converged turbulent en- 
ergies during gravitational collapse. In addition, the turbulent en- 
ergy cascade is expected to generate turbulence even on subgrid- 
scales (SGS), which may backreact on the resolved scales via tur- 
bulent viscosity. In engineering and many other disciplines of com- 
putational fluid dynamics, SGS turbulence models are used to rep- 
resent the effect of unresolved turbulence on resolved scales in so- 
called large eddy simulations (LES). 



3 MAIN RESULTS 

In total, we have performed 9 LES and 9 equivalent ILES for nine 
halos as listed in table [T] for a constant radiation background with 
J 2 i = 10 3 of the H 2 photo-dissociating radiation field for a radiation 
temperature of 10 5 K. The results obtained from our cosmologi- 
cal simulations are presented here in detail. In the early phases of 
the collapse gas falls into the dark matter potential and gets shock- 
heated during the nonlinear evolution phase. Gravitational energy 
of the halo is continuously transferred to kinetic energy of the gas 
and dark matter during the course of virialization. 



3.1 LES Runs 

The properties of the halos at larger scales during the formation 
of the first peak are shown in figure [T] The density profiles fol- 
low an R -2 behavior for all halos according to the expectation 
of an isothermal collapse. The small bump in the density profile 



of one halo shows the formation of an additional clump. The top 
right panel of figure Q] shows the temperature radial profiles. The 
gas in these halos is heated up to their virial temperatures and 
subsequently cools via Lyman alpha radiation. It is found that all 
halos collapse isothermally in the ubiquity of strong UV photo- 
dissociation flux and the formation of H 2 remains suppressed. The 
temperature in the center of the halos is about 8000 K. The dissim- 
ilarities in temperature profiles around radii of 10 9 AU arise from 
the difference in halo masses. The radial profiles of radial velocity 
scaled by the sound speed are depicted in the bottom left panel of 
figure Q] Overall, this ratio is negative which indicates the infall of 
gas to the center of the halo. This behavior is observed for all the ha- 
los. The turbulent Mach number is plotted in the bottom right panel 
of the figure Q] The value of turbulent Mach number is about one 
which indicates that turbulence is trans-sonic in the haloes. We at- 
tribute small variations in the turbulent Mach number to the differ- 
ences in halo masses. The fraction of molecular hydrogen is shown 
in figure |2]f or a few representative cases. The H 2 fraction increases 
at lower densities due to the rise in electron abundance during the 
non-linear phase of the collapse. At intermediate densities, the H 2 
abundance becomes constant as gas cools, recombines and remains 
neutral with a constant temperature around 8000 K. The presence 
of sharp spikes in the H 2 fraction is due to the shocks occurring at 
the central densities due to collisional dissociation as well as due to 
the presence of clumps. In general, the H 2 fraction is lower than the 
universal value (i.e., 10~ 3 ). Therefore, the contribution of H 2 cool- 
ing in the thermal evolution of the halos studied here is negligible. 
The overall halo properties are in accordance with previou s studies 
JAbel et al.ll2002l : IWise et alhoofl iTurket al.ll2009ll2012l) . 

As a result of the strong photodissociating background, the ha- 
los go through an almost isothermal collapse regulated by atomic 
hydrogen line cooling. The central core develops a highly turbulent, 
extended structure. When the highest refinement level is reached, 
the evolution becomes adiabatic, thus stabilizing the gas on the 
finest grid and preventing collapse to smaller scales. The evolu- 
tion on larger scales however continues, leading to turbulent accre- 
tion onto the central clump and occasional fragmentation as seen in 
Fig. [3] In the final stage of the simulation, fragmentation is clearly 
visible in 3 out of 9 halos. Even in some of the other halos, frag- 
mentation occasionally occurs, but the clumps merge again on short 
timescales. Zooming in further, we note that the central clumps are 
surrounded by self- gravitating accretion disks with prominent spi- 
ral arm structures depicted in Fig. |4] The growth of the clumps is 
thus regulated due to accretion via the observed spiral arm struc- 
tures as well as occasional mergers. 

The radial profiles of the physical properties of the halo in the 
central 10 5 AU are shown for four representative cases in figure 
\5\ The mass profile increases as R 3 in the region of constant den- 
sity, is then approximately constant within the disk, thus implying 
a minor contribution of the disk to the total mass, and subsequently 
increases as R according to the density profile. The latter shows R 2 
behavior corresponding to an isothermal collapse, while flattening 
occurs in the central clump, where the collapse is stabilized by an 
adiabatic core. The rotational velocity (V rot = vx ^ where a = Lxr, 
L is the specific angular momentum and r is the radius ) is about 10 
km/s in the center, peaks at radii of 100 AU and then declines as 
it follows the Keplerian rotation. Variations in the peak rotational 
velocity arise due to the difference in the halo mass and spin. The 
radial velocity profile indicates that gas is falling into the center of 
the halos with typical velocity of about 10 km/s. The mean accre- 
tion rate for these halos is thus a few M©/yr. 

Turbulence production and dissipation: We found that the 
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Figure 1. The figure shows the radially binned spherically averaged radial profiles for nine different halos for SGS turbulence. Each line style and color 
represents a halo. The upper left panel of the figure shows the density radial profiles. The temperature radial profiles are depicted in the upper right panel. The 
bottom left panel of the figure shows the radi al velocity profi les scaled by the sound speed. The averaged radial profiles of turbulent mach number (v tur b/C s , 
for the definition of v tur b see equation 4 of JLatif et al.l2013bh ) are depicted in the bottom right panel of the figure. 
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Figure 2. The fraction of molecular hydrogen for two representative cases is shown in the figure. It can be found that this fraction remains lower than universal 
abundance of H2. 



medium inside the halo becomes highly turbulent and clumpy due 
to well resolved turbulence below a scale of ~ 10 5 AU. The tur- 
bulent energy is low in the outskirts of the halo and tends to be 
stronger in the core. The rate of SGS turbulence production is com- 
puted as 






(11) 



while the turbulent dissipation rate pe is defined in equation ([9]). 
The ratio of turbulence production to dissipation for different LES 
is depicted in the left panel of figure [6] In the outer regions, the tur- 
bulent production rate tends to be higher than the dissipation rate. 
This is a consequence of gravitational instabilities exciting turbu- 



lent motions during the collapse. In the adiabatic cores, on the other 
hand, production and dissipation are roughly in equilibrium, which 
indicates developed turbulence. However, local enhancement of 
turbulence production by instabilities can be seen even in the in- 
ner regions of some halos. 

On the basis of a simple Kolmogorov scaling for homoge- 
nous turbulence, one would naively expect the turbulent energy to 
decrease towards small radii. This effect is compensated as turbu- 
lence is produced locally during the gravitational collapse (see also 
lLatif et alj|2013ah . In order to elucidate the contribution of turbu- 
lence, we followed the dynamics of gas compression employing a 
differential equation for the divergence of the velocity field d - V- v 
techmidtetalJ2013h : 
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Figure 3. The figure illustrates the state of simulations at the collapse redshifts for 9 halos. 
for a fixed Jeans resolution of 64 cells. 



Density projections are shown for the central 1000 AU of the halo 



Dd * „ c 

dF =47W ' 



A 



(12) 



Here, ^ = j t +v V, 6 is the overdensity relative to the mean density 
p and A the local support terms against gravitational compression. 
A receives the contributions from thermal pressure, resolved turbu- 
lence and the SGS turbulent pressure. 

Support terms : Thermal pressure support against gravity de- 
fines the critical mass for gravitational instability to become effec- 
tive and is often accounted using the Jeans mass. However, in a 
strict theoretical sense, it is not applicable to substructures in a tur- 
bulent self-gravitating medium. For a more sophisticated analysis, 
we have computed the local support by thermal pressure, turbu- 
lence on numerically resolved scales and SGS turbulent pressure. 
These local support terms are d enoted by A therm , A turb and A S gs ■ 
They are defined as dSchmidt et al .1120131) : 

1 d 2 P 



1 

+ — 



dp dP 
p dxfdxj p 2 dxi dxt 



(13) 



(14) 



where P is pressure, p is the gas density, a> - V X v is the vorticity 
of the fluid and \S |= (2SijSij) l/2 is the rate of strain 



dvj 

dxj 



Su=±l^- + ^ 



dvj_ 
dxj 



Asgs - — 



1 d 2 P SG S 1 dp dPsGS 



p dxidxi p 2 dxi dxi 



(15) 



(16) 



Here P S gs = fp^t is the SGS turbulent pressure. If these source 
terms are positive they provide support against gravity otherwise 
they aid to gas compression. These terms regulate the evolution of 
the divergence of the velocity field given by 

Dd 

- — = AnGp 6 - A. (17) 

The positive and negative contributions of these terms scaled 
by the gravitational compression rate AnGpo5 are plotted in fig- 
ure [6] The plot clearly demonstrates that dynamical compression 
by gravity dominates above 1000 AU because A/AnGpod «: 1. 
The thermal support of the gas becomes dominant and acts against 
gravitational compression at radii smaller than 10 AU, which cor- 
responds to the center of our adiabatic core. In between these 
scales, while the resolved velocity field compresses the gas (neg- 
ative net support by resolved turbulence), this effect is balanced by 
the large positive support from SGS turbulence. In this region, self- 
gravitating disks are able to form due to the additional turbulent 
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Figure 4. Density projections for nine halos (same as figure in the central 300 AU. 



viscosity. Thus, turbulence has both a stabilizing and compressive 
effect, where the former is mainly due to eddies on length scales 
below the resolution limit. The positive support by resolved turbu- 
lence may also include the contribution from rotation. However, it 
is important to note that later is subdominant i.e., the main positive 
support results from SGS model, while resolved turbulence appears 
to have predominantly compressive modes. 

A further analysis of the central regions shows that the disks 
seen in Fig. [4] are rotationally supported and stable against gravi- 
tational instability. The ratio of rotational to Keplerian velocity for 
disks is shown in figure |7] It is found that the disks are partly rota- 
tionally supported as the mean rotational velocity is 0.4 times the 
Keplerian velocity. 

3.2 ILES Runs 

The general properties of the halos for standard hydro runs with- 
out SGS model are identical to those found from the LES (see also 
lLatif et all J2013ah ). All halos collapse isothermally and follow an 
R -2 behavior. The central temperature in the halos before the adia- 
batic evolution is about 8000 K. No differences are observed in the 
thermal properties of the halo. 

We however see strong differences in the morphology and the 



fragmentation of the halos compared to the LES. The state of the 
simulations on a scale of 1000 AU is shown in figure [8] It is found 
that fragmentation is more common in these runs (7 out of 9 halos). 
Nevertheless, we see that a massive object is formed also in the 
core of these halos. The typical masses of these clumps are slightly 
smaller than in the LES in most of the runs due to the presence 
of additional clumps. The masses of the most massive clumps are 
listed in table \T\ The formation of less massive clumps is much 
more common in these runs. We attribute the formation of clumps 
to turbulent compression which locally aids in gravitational com- 
pression and helps in the formation of smaller clumps. Here this 
effect is enhanced as the viscous damping from unresolved tur- 
bulence is missing in these calculations. To quantify the effect in 
more detail, we have compared the contribution of the local support 
terms. They are shown in figure [9]for a representative case. It can 
be noticed that gravitational compression dominates at larger scales 
while the resolved turbulence becomes important in the central 100 
AU, where the latter becomes comparable to the support from ther- 
mal pressure. It may even exceed the thermal support locally. Over- 
all, the support from resolved turbulence is negative which aids in 
gravitational compression. In the very center of the halo, we see an 
adiabatic core which is supported by thermal pressure. 

We also found that in some cases these clumps are produced 
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Figure 5. Radial profiles for the physical properties of the disks centered on the peak density at collapse redshifts. The rotational and radial velocity profiles are 
shown in the top panels. Mass accretion rate (4nR 2 pv ra d) and mass radial profiles are shown in the bottom panel. These profiles are shown for 4 representative 
halos. 
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Figure 6. The left panel of this figure shows the ratio of turbulence production to dissipation for various haloes at their collapse redshifts. Comparison of 
local support terms like thermal pressure, resolved turbulence and SGS turbulence for a representative case is shown in the right panel. They are scaled by the 
gravitational compression. The positive components mean support against gravity while negative components aid the compression. 



at earlier times and merge with the central object. It can be no- 
ticed from the figure that these clumps are very well separated and 
may lead to the formation of binary or multiple systems. Depend- 
ing on the further evolution the clumps might evolve to a massive 
star cluster surrounding a central black hole. We see accretion rates 
of few M Q /yr comparable to LES. Zooming into the central object 
as shown in figure [TO] the clumps are dense and compact. No disks 
are observed in these runs due to the higher turbulent velocities. 
Nevertheless, an object of mass between 500-1000 M© forms in the 



center of every halo which will lead to formation of massive black 
hole or a supermassive star as an intermediate stage. 



4 DISCUSSION 

We present here the highest resolution cosmological large- eddy 
simulations to date which track the evolution of high-density re- 
gions on scales of 0.25 AU beyond the formation of the first peak, 
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Figure 7. The ratio of rotational to Keplerian velocity for the disks shown in figure [4] Each line style and color represents one disk. 




Figure 8. The figure illustrates the state of simulations for various halos (ILES) at their collapse redshifts (at t/f = 4). The density projections are shown for 
the central 1000 AU of the halos. 
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Figure 9. The local support terms for a representative case are shown at collapse redshift (at t/f = 4). The solid lines represent the positive component while 
dashed lines indicate the negative component of the support term. The positive components mean support against gravity while negative components aid the 
compression 



Table 1. Properties of the simulated halos are listed here. 



Model 


Mass 


spin parameter 


Collapse redshift 


Clump Masses 


Clump Masses 




M 


A 


z 


LES (M Q ) 


ILES (Mq) 


A 


8.06 x 10 6 


0.0347468 


12.06 


950 


460 


B 


4.3 x 10 6 


0.0309765 


11.3 


850 


850 


C 


2xl0 7 


0.0178532 


12.6 


800 


611 


D 


1.0 xlO 7 


0.0338661 


12.8 


850 


842 


E 


1.9 xlO 7 


0.0084786 


13.7 


1200 


741 


F 


4.5 x 10 7 


0.0294066 


18.1 


800 


588 


G 


2.3 x 10 7 


0.021782 


15.9 


800 


815 


H 


9.7 x 10 6 


0.0099387 


13.5 


900 


1522 


I 


8.2 xlO 6 


0.0252206 


15.0 


556 


1000 



and study the impact of subgrid-scale turbulence. Our findings 
show that while fragmentation occasionally occurs, it does not pre- 
vent the growth of a central massive object resulting from turbu- 
lent accretion and occasional mergers. The central object reaches 
~ 1000 solar masses within 4 free-fall times. They may further 
grow up to 10 6 solar masses through accretion in about 1 million 
years. The direct collapse model thus provides a viable pathway of 
forming high-mass black holes at early cosmic times. 

The first massive clump is formed in the center of the halo 
as early as z = 18.0, and begins to accrete gas from the ambient 
medium. Due to the additional turbulent viscosity, the morphology 
of the halos becomes more extended compared to the runs without 
SGS model. At the final stage, our simulations harbor central mas- 
sive objects with ~ 1000 M©. These objects are expected to form a 
supermassive star on the Kelvin-Helmholtz timescale, which may 
subsequently collapse to a supermassive black hole dFrver et al.1 
l200ll) . In the presence of efficient accretion, the atmospheres of 
such massive prot ostars remain at low tem peratures, and UV feed- 
back is negligible dHosokawa et al.l2012h . On the other hand, feed- 



back via the accretion luminosity may become relevant at an earlier 
stage. 

Accretion Luminosity : We also computed the heating rate 
from the accretion luminosity for the massive clumps in the center 
of the disk employing the optically thin approximation, as 



^ / L acc \ _j _j 

Tacc = Kp \-^) ergg S 



(18) 



where /c P is the Planckian opacity, R is the distance from the source 
and L acc is the accretion luminosity computed as 



Lace 



GMM 
R* 



(19) 



here M* is the mass of the star, M is the mass a ccretion rate and 
R* is the radius of star obtained from the relation (IHosokawa et al.1 

Im3) 



R* = 2.6 x 10 



'rJ, 



M* 






(20) 



\100M o y 

Our estimates show that for a 500 M clump with a typical size 
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Figure 10. Same as figure[8] density projections for the fields of view of 300 AU, for the ILES runs. 



of R ~ 10 AU, a temperature o f 7500 K, a Planckian opac- 
ity of 0.04 (IMaver & Duschll l2005h and a typical accretion rate 
of 1 M© yr _1 , the heating rate from the accretion luminosity is 
2 x IO -4 erg cm -3 s _1 , which is approximately equal to the atomic 
line cooling, given as 1.6 x IO -4 erg cm -3 s _1 . Once the masses 
of the clumps reach 1000 M Q , the heating rate exceeds the cooling 
rate by a factor of 10. We note, however, that this is easily compen- 
sated by a minor change in the temperature. Specifically, a temper- 
ature rise by 500 K would be sufficient to account for the additional 
heating. However, with the optically thin approximation, we likely 
overestimate the contribution on the scales considered here. In fact, 
with an optical depth r ~ RpK P > 30, the accretion luminosity 
feedback is likely confined to much smaller scales, thus not affect- 
ing the evolution considered here. It can thus be expected that the 
accretion luminosity heating would only cause minor differences 
until this point, and potentially stabilize the central region further, 
i.e., suppress fragmentation even further. 

We further note that th is effect has been examined in further 
detail bv lSmithet al.l (1201 11) for minihalos, finding that the amount 
of fragmentation is at most slightly reduced. We would expect that 
to favor the formation of a central massive object even further. 
However, in the atomic cooling halos considered here, atomic hy- 



drogen is expected to act as an even more efficient thermostat, and 
the relative change of the thermal pressure even less relevant. 

These massive clumps will collapse to similar-mass black 
holes at the end of their lifetime, which can grow to the ob- 
served IO 9 Mq black holes at z ~ 6 - 7 via Eddington accretion 
dShapiroir2005h . The main requirement to obtain high-mass black 
hole seeds is thus a massive primordial halo at z ~ 15, exposed to 
a sufficiently strong radiation field for H 2 dissociati on. The abun- 
dance of such halos was pre viously calculated by iDiikstra et al.1 
d2008l) : lAgarwal et al.l (120121) . showing that these are indeed fre- 
quent enough to explain the number of observed high-redshif t black 
holes. The high-resolution simulations presented here strongly con- 
firm the feasibility of this scenario from a dynamical point of view, 
as a result of turbulent accretion and occasional mergers with addi- 
tional fragments. For a final conclusion on the origin of supermas- 
sive black holes, future studies will need to address the growth of 
these clumps in the presence of feedback. At the later stages, we 
expect star formation and stellar feedback to affect the subsequent 
evolution, potentially giving r ise to a co-evolution between stellar 
bulge and central black hole (JTremaine et al.ll2002c I Johnson et al.l 
l201ll) . The formation of their seeds is in any case feasible, thus 
yielding strong support for direct collapse scenarios. 
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